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Any object on earth has two fundamental properties: it is finite, and it is 
made of atoms. Structural information about an object can be obtained from 
diffraction amplitude measurements that account for either one of these traits. 
Nyquist-sampling of the Fourier amplitudes is sufficient to image single parti- 
cles of finite size at any resolution. Atomic resolution data is routinely used to 
image molecules replicated in a crystal structure. Here we report an algorithm 
that requires neither information, but uses the fact that an image of a natural 
object is compressible. Intended applications include tomographic diffractive 
imaging, crystallography, powder diffraction, small angle x-ray scattering and 
random Fourier amplitude measurements. 
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1. Intro 

In a standard imaging system, light scattered from an 
object forms a diffraction pattern which encodes infor- 
mation about the object Fourier components. A lens 
recombines the scattered rays so that they interfere cor- 
rectly to form an image: it performs an inverse Fourier 
transform of the diffraction pattern to convert Fourier 
(reciprocal) representation of the object into real space 
information. 

At visible wavelengths, aberration-free lenses can pro- 
vide diffraction limited images within a limited depth of 
field. Smaller wavelengths offer much higher resolutions 
by reducing the diffraction limit. X-rays also offer the 
ability to penetrate through thick objects and allow one 
to examine elemental, chemical, or magnetic information 
by exploiting mechanisms such as resonant X-ray scatter- 
ing. Unfortunately, diffraction limited optics are harder 
to come by for X-rays, whose paths are difficult to ma- 
nipulate. For higher resolutions, the optics need to cope 
with light scattered to high angles. Building such optics 
is a difficult technical challenge. Currently, focal widths 
of tens of nanometers are achievable, but lenses capable 
of atomic resolution are so far well beyond reach. 

Diffraction and scattering experiments overcome this 
problem by eliminating the need for any optics. The 
concept is to record the scattering pattern created by 
an object and perform the re-interference normally done 
by a lens numerically instead. Since no optical elements 
are used, aberration free images may be obtained with 
resolutions limited in principle only by the maximum 
momentum transfer which can be achieved. However, 
the resulting image is additionally limited by the com- 
puter's ability to recover the entire image from incom- 
plete Fourier information and loss of phase information. 
The intimate relationship between the phase-front and 




Fig. 1. Simple homometric objectsfl], with the same 
Fourier magnitude, are composed by the convolution be- 
tween two non-centrosymmetric objects. 



the direction of propagation seems to suggest that the 
task of recombining x-rays back at the sample position 
would seem hopeless. It is not so under a surprisingly 
small set of conditions. 

The importance of fine sampling of the diffraction pat- 
tern intensity was recognized at an early stage in x- 
ray crystallography. The observation that Bragg diffrac- 
tion undersamples the diffracted intensity pattern [2j |3] 
was followed by the demonstration that the solutions to 
Nyquist sampled diffraction patterns are almost always 
unique [U |5J El |7] (although one could easily make up 
examples where this is not the case -see Fig. [lj. 

These ideas, along with the development of powerful 
light sources, producing collimated beams of coherent x- 
rays, enabled the development of coherent x-ray diffrac- 
tion microscopy [H] (aka lensless or diffractive imaging). 
This technique aims at imaging, through coherent illumi- 
nation, Fourier amplitude measurements and adequate 
sampling, macroscopic objects such as entire cellular or- 
ganisms [5], or nanoporous aerogel structures [TU]. See 



2 



[11) for a review. 

Diffraction microscopy solves the phase problem using 
increasingly sophisticated algorithms based on the sup- 
port constraint, which assumes adequate sampling. The 
object being imaged is limited within a support region S: 

p(r) = 0, if r f S. (1) 

The sampling conditions required to benefit from the 
support constraint have limited the adoption of projec- 
tion algorithms to other experimental geometries that al- 
low only for sub-Nyquist sampling, most notably Bragg 
sampling from periodic crystalline structures. 

Modern sampling theory however, tells us that that 
Nyquist sampling conditions dictated by the support are 
the worst case scenario for an arbitrary object. In other 
words, Shannon was a pessimist: he did not account for 
the signal structure. Compressive sensing theory tells 
us that the number of measurements are dictated by the 
signal structure rather than it's length. By structured we 
mean that the signal has only a few non-zero coefficients 
when represented in terms of some basis, or can be well 
approximated well by a few non zero coefficients: they 
can be described in terms of a few atoms, a few stars, a 
few wavelet coefficients, or possibly a few protein folds. 

In other words, an object of interest is often sparse or 
concentrated in a small number of non-zero coefficients 
in a well chosen basis, i.e. it can be compressed with no 
or almost no loss of information. The meanings of "well- 
chosen" and "of interest" are slightly circular: A basis is 
well-chosen if it succinctly describes a signal of interest; 
likewise, a signal is of interest if it can be described with 
just a handful of basis elements. 

Since we do not know where these few terms are lo- 
cated, conventional wisdom would indicate that one has 
to first measure the full sample at the desired resolution, 
since overlooking an important component of a signal 
seems almost inevitable if the whole haystack isn't thor- 
oughly searched over. There would seem to be no alter- 
native to processing each signal in its entirety before we 
can compress it and store only the desired information 
(such as the location of the atoms in a molecule). 

But a new theory of "compressive sampling" has shown 
how an image of interest or structured signals generally 
can be reconstructed, exactly, from a surprisingly small 
set of direct measurements. Candes and colleagues [T5] 
have defined a notion of "uniform uncertainty" that guar- 
antees, with arbitrarily high probability, an exact solu- 
tion when the signal is sparse and a good approximation 
when it is compressible or noisy. Their uniform uncer- 
tainty condition is satisfied, among others, by Fourier 
measurements of a sparse real space object. 

The question of whether modern sensing theory is ap- 
plicable to Fourier amplitude measurements was first 
raised by Moravec, Romberg and Baraniuk |13j who pro- 
vide an upper bound sampling condition for the success- 
ful retrieval of a sparse signal autocorrelation, and dis- 
cuss other conjectures with far reaching consequences for 
low-resolution undersampled phase retrieval. Since the 



theory is relatively new and not widely known to the 
phase retrieval community, modern sampling theory is 
briefly reviewed. 

2. A nonlinear sampling theorem 

The notion that a diffraction pattern from a sparse object 
can be reconstructed at sub-Nyquist sampling is not en- 
tirely new. The so called "Direct methods" are routinely 
used for atomic resolution imaging of increasing complex 
molecular structures. Direct methods enforce the condi- 
tion that the resulting molecule is composed of a finite 
number of atoms. The conditions for successful ab-initio 
phase retrieval using these methods are strict: it requires 
(1) atomic resolution and (2) about 5 strong peaks per 
atom. Condition (2) means that the algorithms do not 
scale well with a large number of atoms since the number 
of strong reflections decreases rapidly with the number of 
atoms. 

Here we look for an alternative answer from modern 
sampling theory. 

Suppose that one collects an incomplete set of fre- 
quency samples (amplitude and phase) of a discrete sig- 
nal p(r) of length N. The goal is to reconstruct the 
full signal p given only K samples in the Fourier domain 
where the "visible frequencies" are a subset (of size K) 
of the set of all frequencies {0, • • • , N}. 

At first glance, solving the underdetermined system of 
equations appears hopeless, as it is easy to make up ex- 
amples for which it clearly cannot be done. But suppose 
now that the signal p is compressible, meaning that it 
essentially depends on a number of degrees of freedom 
which is smaller than N. Then in fact, accurate and 
sometimes exact recovery is possible by solving a simple 
convex optimization problem. 

Theorem 2.1. (Candes Romberg and Tao As- 
sume that p is n a -sparse, (e.g. n a atomic charges in 
real space with N resolution elements), and that we are 
given K Fourier coefficients with frequencies selected uni- 
formly at random. Suppose that the number of observa- 
tions obeys K < Cn a logN. Then minimizing l\ recon- 
structs p exactly with overwhelming probability. In par- 
ticular, writing C — 22(6 + I), then the probability of 
success exceeds 1 — N~ s . 

The theorem shows that a simple convex minimization 
will find the exact solution without any knowledge about 
the support, the number of nonzero coordinates of r, their 
locations, and their amplitudes which we assume are all 
completely unknown a priori. 

Following [12 we formulate this more explicitly. The 
algorithm that optimizes the i\ norm: 

min \p\i subject to T p — F ,{k G fi}; (2) 

IpIi = £K*)I (3) 

will find the solution with the correct an- 
swer without knowing the support S (S = 
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{1 if |/?o(a:)| > 0, otherwise}), nor the number of 
non-zero elements 1 1 /Co 1 1 o (llpollo = 

In addition, another remarkable result, is that the con- 
cept of "atomicity" in real space is generalized to spar- 
sity in other basis. We can use, instead of a dictionary 
of point atoms, a dictionary of curves, beams, or a basis 
that describes protein folds using a few terms. Finally, 
since the equation above does not depend strongly on N, 
the number of resolution elements or the basis, we can 
choose redundant basis, with more terms than real space 
resolution elements if it helps to describe the object with 
fewer terms. 

The only requirement is to be able to find the min- 
imum of ||p||i = X)I/°( £C )I- These results have already 
been applied to a number of imaging techniques, but they 
require amplitude and phase of the Fourier coefficients. 
The question that we try to address here is: what are the 
implications to the inversion problem of Fourier magni- 
tude only recordings? 

The following theorem gives an upper bound: 

Theorem 2.2. (Moravec, Romberg and Baraniuk f!3f): 
Assume that p is n a -sparse, then it can be recovered 
exactly from a reduced number of random Fourier mag- 
nitude measurements K > Cn\ log(4M/n^). 

The theorem is based on the fact that the autocorre- 
lation of an n a -sparse object is at most n\ sparse. The 
authors also point out that n\ sparse is worst case sce- 
nario. If the object is connected, then the sparsity of 
the autocorrelation grows linearly, instead of quadrati- 
cally, with the object support to twice the object size. 
In summary, this theorem provides an upper bound for 
the exact recovery of the object autocorrelation for an 
arbitrary but finite object, as well as from sub-sampled 
data for an object made of a few atoms. 

The important question then is: does a typical geom- 
etry used in x-ray diffraction satisfy the conditions re- 
quired to enable the full recovery of the object's auto- 
correlation? One of the keys to compressive sensing is 
the role played by randomness in the data acquisition 
(from the Uniform Uncertainty Principlc[12 ). Does a 
Bragg geometry satisfy these conditions? In this section 
we explore the possibility that various sampling geome- 
tries satisfy the condition for compressive recovery of the 
autocorrelation . 

We simulate an object of n a atoms, and attempt to 
recover the full autocorrelation from a subset of Fourier 
amplitude measurements. The optimization problem can 
be stated as follows: 

min \p\i subject to \Ta\ 2 = Ik (4) 

where a is the autocorrelation of the object we are trying 
to reconstruct. 

The problem setting was written within the SPARCO 
toolbox[14j. With the addition of radial averaging, the 
SPARCO toolbox provided the operators to simulate a 
general scattering experiment: a Fourier transform, and 



a Fourier mask. The SPGL1 software [TS] used here was 
able to converge quickly to a root. Various sampling con- 
ditions are explored in Fig. [2] We can see that random 
Fourier sampling and limited angle tomography satisfy 
the conditions for exact recovery, while Bragg sampling 
(sampling every other 2 Fourier components in each di- 
mension) does not. However we are able to recover the 
aliased autocorrelation from a severely reduced number 
of measured Bragg reflections or even from radially aver- 
aged powder data (using an oblique unit cell which causes 
a reduced number of peak overlaps, in this case an aver- 
age of 10 peak overlap over the same radial shell). 

3. Compressive phase retrieval algorithms 

Although results described above are encouraging, it is 
obvious that we can do better than that: we have not 
utilized the notion that the function that we have been 
trying to recover is itself the autocorrelation of an even 
sparser object. We therefore look for an algorithm that 
recovers a sparse object from subsampled Fourier ampli- 
tude measurements. 

Though problems of this (non-convex) nature are dif- 
ficult because only exhaustive algorithms can guarantee 
convergence, a variety of heuristics perform well in prac- 
tice. We follow the notation described elsewhere [28] for 
projection algorithms. 

A. Projection algorithms 

The aim of projection algorithms is to find a signal that 
lies in the intersection of two constraint sets. The first 
set is that of the measured Fourier magnitudes. To com- 
pute the projector corresponding to the Fourier magni- 
tude constraint, one first needs to propagate p(r) to the 
data space by a Fourier transform J 7 , and then replace 
estimated amplitudes \p\ (p(k) — Tp(r)) with the mea- 
sured ones, and propagate back to real space. Formally, 
we incorporate the forward T and inverse Fourier 
transform in the operator defined in Fourier space P m : 

Pm = J~ PmJ~ i (5) 

and enforce the condition that the image Fourier mag- 
nitudes are equal to the measured ones. Using these 
transforms one simplifies the calculation of the projec- 
tion which becomes an element-wise operation on each 
recovered Fourier component: 

P m ~p(k) = Vm^J-y (6) 

This projection requires each Fourier measurement to be 
sampled individually. In powder diffraction, partial co- 
herent illumination, broadband illumination, the projec- 
tion operator needs to be generalized. If one defines the 
averaging operator A, the generalization is as follows |30j : 
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(a)Sparse object (b)Diffr. Patt. 





(c) Autocorrelation 




(d)Limited Angle (e)Masked FT of (b) 
Tomoe. 



(f)Recovered 
autocorr. 




(g)Random Fourier (h)Masked FT of (b) 
Measurements 



(i)Recovered 
autocorr. 



where Ua is the corresponding value of the intensity 
which \p(k)\ contributes to. Partial overlap can be ana- 
lyzed using the Richardson-Lucy deconvolution [THl E] 
of the intensities under square root (see also 18 1). The 
second projector, used in diffractive imaging, is that of 
the support: P s p(r) ={p{r) if r e S; otherwise} that 
acts element-by-element to the real space basis. 

While early experiments relied on low resolution imag- 
ing to determine the object support, the development 
of automated support refinement techniques |19j have 
enabled diffractive imaging to solve structures indepen- 
dently. Note that these Shrinkwrap support finding- 
techniques rely on the most compact (and sparse) object 
that satisfy the measurements. 

B. Sparsifying algorithms 

Since we do not know the support of our signal, we re- 
place the projection onto the support set P s with a va- 
riety of operators O known to promote sparsity in phase 
retrieval problems: thresholding, Sayre' squaring opera- 
tor [2U], and a soft thresholding operator that promotes 
sparsity through l\ optimization [34 (Moravec et al. sug- 
gest using the £1 norm as a constraint, but this is rarely 
known) : 




(j)Bragg Sampling 



(k) Masked FT of 
(b) 



(1) Recovered AC 




(m)Random Bragg (n)Masked FT of (b) (o)Recovered AC 
sampling 




(p)Radial averaging 
sampled with 
oblique lattice 



(q)FT of radial avg. 
of (b) 



(r) Recovered AC 



Fig. 2. Recovering a sparse autocorrelation. Top 

row: (a) Original object, (b) Fourier magnitude (c), 
autocorrelation. In the successive rows, Fourier mask 
(left), autocorrelation from masked Fourier amplitudes 
(center) , recovered autocorrelation using l\ minimization 
(right). Fourier masks: (d) limited angles (g) random (j) 
Bragg (m) random Bragg (p) radial average -oblique unit 
cell. 



O t x 
2 x 



x if |x| > r, 
otherwise. 



x ~ T R if M > T ' 
otherwise. 



(8) 

(9) 
(10) 



With slight abuse of notation, we use the same symbol 
for a projector, and define a reflector operator: 

Pr,2,il = Tj 2,ei', Rt : 2,£1 = 20 Tj 2,« _ 1 5 

The first most obvious algorithms are simple alternat- 
ing type: 

= P m O (T/h 2)P (n) 

We note that the squaring method obtained by alter- 
nating the operators P m 02X is equivalent to the "tan- 
gent formula" of Direct methods, where the squaring op- 
eration equation in real space is performed directly in 
Fourier space through an autocorrelation. Remarkable 
improvements in the range of convergence has been ob- 
tained by increasing the step produced by these algo- 
rithms in real space by 2, using the reflector operator R. 
The charge flipping algorithm replaces the thresholding 
operation by with an operator that moves twice as far: 



P 



(n+l) — pvt. 



(t,£L,2)P 



(n) 



In appendix we describe why it is superior to simple al- 
ternating projections. 
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Table 1. Summary of various algorithms 



Algorithm 


Iteration p (n+1) = 


ER 


7-1 i~j (n) 
PsPmP ' 


SF 


i~> 7~» ( n) 

R s P m p y ' 


HIO 


fp ra p< n >(r) reS 

1/7" n n \ { n\ / \ n 

[\I — (iP m )p ' (r) rfS 


DM 


17 + BP, Ml +t„) P™ — -y.II 

- /3P m [(l + 7m )P. - 7m 7]}p (n) 


ASR 


|[P s P m +/]p (n) 


HPR 


§[«, (JJm + (/3-l)Pm) 

+1 + (1 - /3)P m ]p (n) 


RAAR 


[|/3 (R s R m +!) + (!- /3)P m ] P (n) 



However we arc interested in an algorithm that is more 
general than this for reasons that will become apparent in 
the following section. In particular we need an algorithm 
that is robust against the relaxation of the positivity of 
the object, and the change of basis used to describe the 
object. 

The algorithms tested include: HIO [22], SF(23], DM 
[23] , HPR[25] and RAAR [25] (see [28] for a review). We 
tested algorithms based on these ideas for increasingly 
complex phase retrieval problems. A unit cell in a peri- 
odic system was filled with a limited number of atoms. 
First an algorithm has to be stable around the solution. 
If perturbed from the solution, it should go back or at 
least not diverge much from it. Perturbations tested in- 
cluded: distributed noise, "salt and pepper" noise and 
a single large extra charge added to the structure. Sec- 
ondly, it should converge to the solution starting from an 
arbitrary set of phases for a large number of atoms n a . 
For small n a all but a handful of algorithms work, as n a 
increases, it takes longer to converge, and one by one, 
various algorithms stop working. The best algorithm is 
the one that works with highest ratio between number of 
atoms and number of measurements. 

Through these numerical tests, the following outper- 
formed all others (fig. 3): 

51 = \P m P ~ p\> Ti, 

5 2 = P mP >0 

= (SikS 2 ) (P m p {n) ) + (1 ~ Si) 

- (3PmP {n) ; (11) 

where > is intended as a relational operator, Si ; 2 are 
binary masks. If positivity cannot be enforced, then S2 — 
1. From here on, this algorithm will be referred to as 
Espresso, in honor of the compressed coffee. 

4. Beyond atoms 

The algorithms described perform remarkably well in 
case an object being imaged is sparse in real space, 
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Fig. 3. Recovering a sparse object in real space. 

Comparison of various algorithms. Top: test unit cell, 
64x64 resolution elements, 160 atoms of values ranging 
from 1 to 3 (arbitrary units). Bottom: success rate vs 
number of iterations, using the optimal threshold level 
(r = .5 for charge flipping, .4 for Espresso); and (3 = .5. 



namely that it is composed of a limited number of atoms. 
Object made of more than a few atoms, such as large 
macromolccules or biological cells, require a different ap- 
proach. 

Natural objects are often characterized by the fact that 
they are compressible in some basis. Everyone is familiar 
with the fact that images of millions of pixels can be 
saved with nearly indistinguishable accuracy at a small 
fraction of the initial image size. In other words, natural 
objects can often be accurately described in terms of only 
a few non-zero coefficients in some basis. 

Moravec et al.[13] conjecture that compressive meth- 
ods for phase retrieval could be applicable in a basis other 
than the real space. Here we set out to test this con- 
jecture for the crystallographic phase problems with the 
Espresso algorithm. Will all basis work? not if they are 
too localized in Fourier space, as preliminary tests with 
curvelets [3T] seem to suggest. 

The problem in this test is to recover a Schepp-Logan 
phantom test object. As in many tomographic settings, 
we seek to sparsify the the gradient of the object. 

The discrete gradient is obtained by simple matrix op- 
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erations. One defines a two-dimensional matrix D x of 
the same size as p, with the first two elements equal 
to -1 and +1 respectively, and the rest of the elements 
equal to 0. The discrete derivative is obtained by a 
simple convolution between the object and this matrix: 
d x p(r) — D x * p(r). Convolution becomes a product 
in Fourier space, where we define the discrete Fourier 
transform D x — TD X : d x p — J-^ 1 {D x p). Similar argu- 
ments apply for the other direction. We define the matrix 
D y = D x as the transpose of D x . 

The inverse operation used to recover the object from 
the compressed gradient values (d x p, d y p) was the follow- 
ing peudoinverse: 



p = T- 



D\Td yP 



btb x 



DyDy 



(12) 



with e a regularization term. Once we have obtained the 
gradient, [d x ,d y ]p, we apply the Espresso algorithm in 
this space, trying to compress the gradient. The norm 
used was y/\d x p\ 2 ■ 



\dyPf 



Si 

s 2 

d xP n+1 

dyP n+1 



\d x P m p - p\ 2 + \d y P m p - p\ 2 > n, 



although the two gradient components could be treated 
separately, each with their own supports Si t 2- 

Examples are shown in Fig. 4, where different masks 
are applied to the Fourier data: Full sampling, missing 
central region due to the beamstop, limited angles of to- 
mographic datasets, Bragg sampling, and radially aver- 
aged powder data. Notably, the periodic crystal symme- 
try is "discovered" by sampling the diffraction pattern 
at the Bragg condition. In the case of powder diffrac- 
tion, with an average of 20 overlaps per sampled point, 
the algorithm is still able to recover the outline of the 
object, or the dominant terms contributing to the gra- 
dient. The unit cell used for calculating the radial aver- 
age was oblique, to reduce some peak overlaps. In each 
case the threshold is defined dynamically: r — 2a, with 
a = rms(\p- P m p\). 

5. Conclusions 

The revolutionary findings of compressive sensing have 
shown how severely incomplete data can be used to re- 
cover both sparse and compressible signals from Fourier 
amplitude and phase measurements. This paper dis- 
cusses its implications for the phase retrieval problem, 
namely that one does not need to assume a finite signal 
bandwidth, nor atomic resolution. A new algorithm for 
compressive phase retrieval has been proposed. It outper- 
forms existing algorithms for atomic resolution data, and 
can be applied at arbitrary resolution to crystal diffrac- 
tion, provided that the sample is compressible in some 
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Fig. 4. Recovering an object without atomic res- 
olution. Shepp-Logan phantom is embedded into a 
128x128 image. Left: Fourier mask, Right: reconstruc- 
tion. (aWb): Full sampling, (cWd) Beam stop, (e)-(f) 
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basis. Optimization of the parameter j3 on the fly at 
each iteration |32j . promises to further improve this al- 
gorithm and will be subject of future work. The search 
for the optimal sparsifying basis and improved compres- 
sive algorithms promises to revolutionize low resolution 
phasing in X-ray diffraction. 
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Appendix A: Charge nipping 

Take a known object with charge density distribution 
p(r), with Fourier amplitude \p(k)\ = ^I{k). We first 
look at perturbations from the solution. If we add a ran- 
dom noise distribution a with root mean square (r.m.s.) 
(j < min(p : p > 0), a simple threshold will be able to iso- 
late the signal p > a from the noise. As the noise grows, 
an increasing large percentage of noise will survive the 
thresholding operation and the thresholding algorithm 
will stagnate at first local minimum it encounters. At 
low perturbation levels, a simple thresholding operation 
is the obvious and effective answer. Flipping components 
below a threshold, rather than setting them to zero would 
just flip sign within a first order term, and at best remove 
second order 0(a 2 /||p||) contributions to the noise a. We 
consider "salt and pepper" noise a = ^i=i a t^( x ~ x i) 
with a few s terms, such as if we displaced a small per- 
centage of the peaks from the correct answer) . For sim- 
plicity, let us consider just a single charge a in x : 



using Taylor series 



l\+x 



1 



iX 



p' = p + a5(x — x ); 



(Al) 



We enforce the measured data by applying the Fourier 
magnitude projector to the current guess p' = T p' = 
p+ae lkx °: 

The projection is performed by Fourier transforming 
p' = T ' p' ', and enforcing the Fourier magnitude: 



PrnP' 




a ikx 



= P 
= P 



hae ikxo 



\ae lkx ° 



1 2|/5| 2a ' 



1 PP 



ikxo 
ikxo i /~\f Jl 



1 P n * r -ikx 

2 |p| 2 " e 



+ 0(a 2 )) 



2p! ut + 0{a ) 
la*F(x ) + O(a 2 ) 



2 



The term has a rapidly oscil- 



with F(xq) = 

lating phase, and diffuses half of the charge a everywhere 
in the unit cell through the function / = T~ l F . At first 
order in a one obtains: 

P m p =P+ \aS{x - as ) ~ la*f(x ) + 0{a 2 ) (A2) 

with f(xo) being distributed over the unit cell. Since the 
charge is reduced only by half, we could try moving twice 
as far, instead of moving from rho to P m rho by a step 
/ — P m , using a reflector R = 2P — I to: 



Rmp' = P - a*f(x ) + 0(a 2 ) 



(A3) 



ore 



-ikx 



Ip! 2 



Suppose the threshold only picks up the term 
|a*/(xo), then the flipping operation converges to the 
solution in one step (at first order): 

P m [p+ \aS{x - z ) + \a*f{x Q )\ = p + 0(a 2 ) (A4) 

Although this result describes only local convergence, it 
shows how a few wrong peaks can be recovered easily. 
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